Computing the Convex Hull
Introduction
Imagine that you have a set of points and you want to draw a polygon so that all of these points are either on the polygon, or inside of it. If you want to make that polygon convex - i.e., if you take a look at it from the outside, no outer angle is lower than 180 degrees -, then you can do that by computing the convex hull of this set of points.
A convex combination is a linear combination of the kind
t * A + u * B + v * C + ...,
where the sum of all coefficients t + u + v + ... = 1. (Without this condition, it would be called an affine combination.) The convex hull of a set of points is the set of all convex combinations of these points with all possible values of the coefficients.
If you have just two distinct points A and B, the convex hull describes the line segment between these two points. If you have three distinct points A, B and C, the convex hull describes a triangle with A, B and C being the corners (unless all three points are on one line, then it's just a line again). With four or more distinct points, it's more complicated. It may be a quadrangle as well as just a triangle or even a line if all of the four points are on the same line again. You will get a triangle if one of the points is inside the triangle formed by the other three points. That's the special property of convex hulls: You won't get polygons that have concave outer angles.
So, basically, the convex hull of a set of points S is the set of all points either on or inside a convex polygon that has all points of S either on its boundary lines or inside. Being able to compute such polygons may be useful for demo effects, so how can it be done? The answer is that there are several algorithms for computing the convex hull. One of them (not the fastest one, but pretty fast considering how easy it is to explain how it works) is the monotone chain algorithm, which has a run time complexity of O(n * log n).
Monotone chain algorithm
To use this algorithm, all points must be ordered by their coordinates. Start with one point that definitely is on the outer line delimiting the convex hull (e.g. the leftmost one). Then add the next two points to the solution set and, in case you started with the leftmost point (if you started with another, you will have to adapt this procedure), check whether they make a counter-clockwise rotation. If so, discard the point in the middle from the solution set - it is inside the polygon. Continue adding the next point and making this check until all points have been processed this way.
By this, you have computed one half of the outer boundaries of the convex hull. Now you have to compute the other half. From the most recently added point, walk back (e.g. to the left) and do everything as previously described, but check for the other type of rotation. As soon as you reach the starting point, you have computed the convex hull.
How to determine whether it's a clockwise or counter-clockwise rotation? Imagine you have three points which are called from left to right A, B and C and you are computing the upper part of the hull. Then you can compare the steepnesses of the lines AB and BC, i.e. you compare (by - ay) / (bx - ax) with (cy - by) / (cx - bx). If the former line is steeper, this means that there is a clockwise rotation, so it is okay to keep the point B inside the set of solutions.
Instead of checking whether
(by - ay) / (bx - ax) > (cy - by) / (cx - bx),
you can, of course, check whether
(by - ay) * (cx - bx) > (cy - by) * (bx - ax).
This will run faster on a CPU on which multiplication is less expensive than division.
Code
/*
Computational Geometry - Convex Hull
Implemented by Claus-Dieter Volko
*/
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <GL/glut.h>
#define NUMPOINTS 20
struct point
{
int x, y;
} points [NUMPOINTS];
int hull [2 * NUMPOINTS + 1];
int *lower;
int *ptr;
int pointCmp (const void *p1, const void *p2)
{
int temp = ((struct point *)p1)->x - ((struct point *)p2)->x;
if (temp < 0) return -1;
if (temp > 0) return 1;
temp = ((struct point *)p1)->y - ((struct point *)p2)->y;
if (temp < 0) return -1;
if (temp > 0) return 1;
return 0;
}
void displayFunc ()
{
int i;
gluOrtho2D (0, 100, 100, 0);
glClear (GL_COLOR_BUFFER_BIT);
glColor3f (1.0, 1.0, 0.0);
for (i = 0; i < NUMPOINTS; i++)
{
glBegin (GL_POINTS);
glVertex2f (points [i].x, points [i].y);
glEnd ();
}
glColor3f (1.0, 1.0, 1.0);
for (i = 0; hull [i] != 0 || i == 0;)
{
glBegin (GL_LINES);
glVertex2f (points [hull [i]].x, points [hull [i]].y);
i++;
glVertex2f (points [hull [i]].x, points [hull [i]].y);
glEnd ();
}
glFlush ();
}
void main ()
{
int i;
int *tptr;
for (i = 0; i < NUMPOINTS; i++)
{
points [i].x = rand () % 99 + 1;
points [i].y = rand () % 99 + 1;
}
qsort (points, NUMPOINTS, sizeof (struct point), pointCmp);
printf ("Points:\n");
for (i = 0; i < NUMPOINTS; i++)
printf ("A%d: (%d, %d)\n", i, points [i].x, points [i].y);
getch ();
ptr = hull;
*ptr = 0;
for (i = 1; i < NUMPOINTS; i++, ptr++)
{
*(ptr + 1) = i;
while (ptr > hull &&
(points [*(ptr - 1)].y - points [*(ptr + 1)].y)
* (points [*ptr].x - points [*(ptr - 1)].x) >=
(points [*(ptr - 1)].y - points [*ptr].y)
* (points [*(ptr + 1)].x - points [*(ptr - 1)].x))
{
*ptr = *(ptr + 1);
ptr--;
}
}
printf ("\nUpper hull:\n");
for (tptr = hull; tptr < ptr; tptr++)
printf ("A%d, ", *tptr);
printf ("A%d\n", *tptr);
getch ();
lower = ptr;
for (i = NUMPOINTS - 2; i >= 0; i--, ptr++)
{
*(ptr + 1) = i;
while (ptr > lower &&
(points [*(ptr + 1)].y - points [*(ptr - 1)].y)
* (points [*(ptr - 1)].x - points [*ptr].x) >=
(points [*ptr].y - points [*(ptr - 1)].y)
* (points [*(ptr - 1)].x - points [*(ptr + 1)].x))
{
*ptr = *(ptr + 1);
ptr--;
}
}
printf ("\nLower hull:\n");
for (tptr = lower; tptr < ptr; tptr++)
printf ("A%d, ", *tptr);
printf ("A%d\n", *tptr);
getch ();
glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
glutInitWindowPosition (50, 100);
glutInitWindowSize (400, 300);
glutCreateWindow ("Computational Geometry - Convex Hull");
glClearColor (0.0, 0.0, 1.0, 0.0);
glMatrixMode (GL_PROJECTION);
glutDisplayFunc (displayFunc);
glutMainLoop();
}